/*
 * case.h
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#ifndef CASEPAM_VTG_H_
#define CASEPAM_VTG_H_

// BASIC LIBLARIES
#include<cstdlib>
#include<cstdio>
#include<iostream>
#include<math.h>
#include<time.h>
// MY LIBRARIES
#include "mros-v42.h"
#include "rootGSL.h"

namespace case_PAM_Vtg
{

const double LANDE_FACTOR_G = 2.;
const unsigned N_PAR = 13; // including L
const unsigned N_PAR_INT = 1;
const unsigned N_VAR = 10;
const char OUTPUT_FILE[]= "test.dat";
const char OUTPUT_FILE_GNUPLOT[] = "test.gnu";

class Par : public Parameters {
	private:
		int NtabD; // counter for tabD double table
		int NtabI; // counter for tabI integer table

		double *tabD;
		int *tabI;

		enum labelD {_n, _t, _tPrim, _U, _V, _h, _epsilonF, _beta,  _tF, _gFg, _Vint, _alpha};
		enum labelI {_L};

	public:

		Par() : NtabD(-1), NtabI(-1) {tabD = NULL; tabI = NULL;}
		Par(unsigned n) : NtabD((int)n - N_PAR_INT), NtabI(N_PAR_INT), tabD(new double[NtabD]), tabI(new int[NtabI]) {
			zeroAll();
		}
		~Par(){
			if(NtabD != -1) { delete[] tabD; NtabI = -1; }
			if(NtabI != -1) { delete[] tabI; NtabD = -1; }
		}
		Par(Par& other) : NtabD(other.NtabD), NtabI(other.NtabI) {
			if(NtabD > 0) {
				tabD = new double[NtabD];
				for(int i=0; i<NtabD; i++) { tabD[i] = other.tabD[i]; }
			}
			if(NtabI > 0) {
				tabI = new int[NtabI];
				for(int i=0; i<NtabI; i++) { tabI[i] = other.tabI[i]; }
			}
		}

		void setNoOfDoubles(unsigned n) { if(NtabD > 0) { releaseDoubles(); } NtabD = (int)n; if(n > 0 ) { tabD = new double[n]; } }
		void setNoOfInts(unsigned n) { if(NtabI != -1) { releaseInts(); } NtabI = (int)n; if(n > 0 ) { tabI = new int[n]; } }
		void releaseDoubles() { if(NtabD != -1) { delete[] tabD; NtabD = -1; } }
		void releaseInts() { if(NtabI != -1) { delete[] tabI; NtabI = -1; } }

		Par &operator=(const Par &rhs) {
			releaseDoubles();
			NtabD = rhs.getNtabD();
			releaseInts();
			NtabI = rhs.getNtabI();
			if(NtabD > 0) {
				tabD = new double[NtabD];
				for(int i=0; i<NtabD; i++) { tabD[i] = rhs.tabD[i]; }
			}
			if(NtabI > 0) {
				tabI = new int[NtabI];
				for(int i=0; i<NtabI; i++) { tabI[i] = rhs.tabI[i]; }
			}
			return *this;
		}

		double n() const { return tabD[_n]; };
		double t() const { return tabD[_t]; };
		double tPrim() const { return tabD[_tPrim]; };
		double U() const { return tabD[_U]; };
		double V() const { return tabD[_V]; };
		double h() const { return tabD[_h]; };
		double epsilonF() const { return tabD[_epsilonF]; };
		double beta() const { return tabD[_beta]; };
		double tF() const { return tabD[_tF]; };
		double gFg() const { return tabD[_gFg]; };
		double Vint() const { return tabD[_Vint]; };
		double alpha() const { return tabD[_alpha]; };
		int L() const { return tabI[_L]; };

		int getNtabD() const { return NtabD; }
		int getNtabI() const { return NtabI; }

		void set_n(double d) { tabD[_n]=d; };
		void set_t(double d) { tabD[_t]=d; };
		void set_tPrim(double d) { tabD[_tPrim]=d; };
		void set_U(double d) { tabD[_U]=d; };
		void set_V(double d) { tabD[_V]=d; };
		void set_h(double d) { tabD[_h]=d; };
		void set_epsilonF(double d) { tabD[_epsilonF]=d; };
		void set_beta(double d) { tabD[_beta]=d; };
		void set_tF(double d) { tabD[_tF]=d; };
		void set_gFg(double d) { tabD[_gFg]=d; };
		void set_Vint(double d) { tabD[_Vint]=d; };
		void set_alpha(double d) { tabD[_alpha]=d; };
		void set_L(int n) { tabI[_L]=n; };

		void print(FILE *file=stdout) const;
		void print_verboten(FILE *file=stdout) const;
		void print_label(FILE *file=stdout) const;

		void zeroAll() {
			for(int i=0; i<NtabD; i++) { tabD[i] = 0.; }
			for(int i=0; i<NtabI; i++) { tabI[i] = 0; }
		}
};

class Var : public Variables {

	private:
		int NtabD; // counter for tabD double table

		double *tabD;

		enum labelD {_lnF, _lmF, _nF, _mu, _mF, _dF, _chiup, _chidown, _lchiup, _lchidown};

	public:

		Var() : NtabD(-1) {tabD = NULL;}
		Var(unsigned n) : NtabD((int)n) { tabD = new double[n]; zeroAll(); }
		~Var() { if(NtabD != -1) { delete[] tabD; NtabD = -1; } }
		Var(Var& other) : NtabD(other.NtabD) {
			if(NtabD > 0) {
				tabD = new double[NtabD];
				for(int i=0; i<NtabD; i++) { tabD[i] = other.tabD[i]; }
			}
		}
		double &operator[](int i) {
			if(i < NtabD && i >= 0) return tabD[i];
			else exit(EXIT_FAILURE);
		}
		const double &operator[](int i) const {
			if(i < NtabD && i >= 0) return tabD[i];
			else exit(EXIT_FAILURE);
		}
		Var &operator=(const Var &rhs) {
			releaseDoubles();
			NtabD = rhs.getNtabD();
			if(NtabD > 0) {
				tabD = new double[NtabD];
				for(int i=0; i<NtabD; i++) { tabD[i] = rhs.tabD[i]; }
			}
			return *this;
		  }

		void setNoOfDoubles(unsigned n) { NtabD = (int)n; if(n > 0 ) { tabD = new double[n]; } }
		void releaseDoubles() { if(NtabD != -1) { delete[] tabD; NtabD = -1; } }

		double lnF() const { return tabD[_lnF]; }
		double lmF() const { return tabD[_lmF]; }
		double nF() const { return tabD[_nF]; }
		double mu() const { return tabD[_mu]; }
		double mF() const { return tabD[_mF]; }
		double dF() const { return tabD[_dF]; }
		double chiup() const { return tabD[_chiup]; }
		double chidown() const { return tabD[_chidown]; }
		double lchiup() const { return tabD[_lchiup]; }
		double lchidown() const { return tabD[_lchidown]; }
		//int getNtabD() { return NtabD; }
		int getNtabD() const { return NtabD; }

		void set_lnF(double d) { tabD[_lnF] = d; }
		void set_lmF(double d) { tabD[_lmF] = d; }
		void set_nF(double d) { tabD[_nF] = d; }
		void set_mu(double d) { tabD[_mu] = d; }
		void set_mF(double d) { tabD[_mF] = d; }
		void set_dF(double d) { tabD[_dF] = d; }
		void set_chiup(double d) { tabD[_chiup] = d; }
		void set_chidown(double d) { tabD[_chidown] = d; }
		void set_lchiup(double d) { tabD[_lchiup] = d; }
		void set_lchidown(double d) { tabD[_lchidown] = d; }

		void set(Variables *V) {
			for(int i=0; i<NtabD; i++) { tabD[i] = (*V)[i]; }
		}

		void shake(double d=1e-8) {
			for(int i=0; i<NtabD; i++) { tabD[i] += (2.*frand()-1.)*d; }
		}

		void print(FILE *file=stdout) const;
		void print_verboten(FILE *file=stdout) const;
		void print_label(FILE *file=stdout) const;

		void zeroAll() {
			for(int i=0; i<NtabD; i++) { tabD[i] = 0.; }
		}
};

typedef struct Global {
	double Energy;
	double magnetization;

	double qup;
	double qdown;
	double q2up;
	double q2down;
} global;

typedef struct Debug {
	double fq;
	double fdqdnF, fdqdmF, fdqddF;
	double fE;
	double fdEdlnF, fdEdlmC, fdEdlmF, fdEdq;
} debug;

class Res : public Result {
	private:
		int n;
		Variables *V;

	public:
		Res(int n) : n(n) { if(n>0) { V = new Var(n); } else { exit(EXIT_FAILURE); } }
		~Res() { delete V; }
		global glob;

		void set(gsl_multiroot_fsolver *s);

		void print(FILE *file=stdout) const;
		void print_verboten(FILE *file=stdout) const;
		void print_label(FILE *file=stdout) const;

		Var* get_V(){ return (Var*)V; }
};

// FUNKCTION (only declaration)
int function(const gsl_vector * x, void *params, gsl_vector * f);

//////////////////////////////////////////////

void example();
void example_inf();
void run_1D();
void run_2D();
void run_Belitz();
void run_Belitz_automatic();
void run_Belitz_2ways();
int automatic_Belitz_2ways(double startg, double starth, double startT, double startHM, double startPM, int startL);
void run_Belitz_automatic2D();
void run_Belitz_2ways_v2();
int automatic_Belitz_2ways_v2(double startg, double starth, double startT, double startHM, double startPM, int startL);


void run_2D_VandVint();
void run_2D_hVint();
void run_2D_hV();
void run_2D_realistic();

//////////////////////////////////////////////

void automatic_run_n(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_U(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_V(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_epsilonF(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_h(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_beta(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_T(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_L(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_tF(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_Vint(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_run_gFg(Par &inpar, Var &invar, double start, double stop, double step);

void automatic_run_VtoP(Par &inpar, Var &invar, double start, double stop, double step, double A, double B);

void automatic_2Drun_nV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_Vn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_hV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_hVtF(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double tFstart, double tFstop, double tFvalue);
void automatic_2Drun_efV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_Vef(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_VintV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_VVint(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);

void automatic_2Drun_tFV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_VtoPn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B);
void automatic_2Drun_nVtoP(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B);
void automatic_2Drun_hVtoP(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B);

void automatic_run_p(Par &inpar, Var &invar, double start, double stop, double step);
void automatic_2Drun_ptF(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2);
void automatic_2Drun_hp(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double tFstart, double tFstop, double tFvalue);
void automatic_2Drun_ph(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double tFstart, double tFstop, double tFvalue);


const double adapt_min_step = 0.00001;
const double adapt_max_step = 0.1;
const double adapt_step_reduction = 2.; // 1.3^3 \approx 2, 1.2^4 \approx 2.
const double adapt_min_m_difference = 0.0001;
const double adapt_max_m_difference = 0.002;
const double adapt_min_difference = 0.001; // if smaller difference between variables, I want to have bigger step
const double adapt_max_difference = 0.02; // if bigger difference between variables, I want to have smaller step
const double coeff_T = 1; // 388.66; // 1 just for inversion of the temperature
const double coeff_H = 577.925;
void adaptic_2D_PMHM_nVtoP(Par &par, Var &initV_PM, Var &initV_HM, double start_n, double stop_n, double step_n, double start_V, double stop_V, double step_V, double A, double B);

} // end of namespace case_PAM_Vtg

#endif /* CASEPAM_VTG_H_ */
